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We consider the density response of a spin-imbalanced ultracold Fermi gas in an optical lattice 
in the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state. We calculate the collective mode spectrum 
of the system in the generalised random phase approximation and find that though the collective 
modes are damped even at zero tempererature, the damping is weak enough to have well-defined 
collective modes. We calculate the speed of sound in the gas and show that it is anisotropic due to 
the anisotropy of the FFLO pairing, which implies an experimental signature for the FFLO state. 
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I. INTRODUCTION 

Ultracold Fermi gases are dilute systems of fermionic 
atoms cooled down to temperatures where quantum 
statistics dominates the physics. The unprecedented ex- 
perimental possibilities of controlling and tuning the ul- 
tracold Fermi gas systems have made them an extremely 
succesful tool for simulating a broad range condensed 
matter phenomena [1-3]. To highlight the topic area of 
this paper, the ability to control the number of each atom 
species forming the ultracold gas has enabled the exper- 
imental study of spin-population imbalanced fermionic 
superfluidity 

One candidate for the theoretical description of imbal- 
anced fermionic superfluids is the Fulde-Ferrell-Larkin- 
Ovchinnikov (FFLO) state pj]-[l3| originally derived for 
superconductors in a strong magnetic field. In the FFLO 
state the pairing correlations which give rise to superflu- 
idity occur with a finite center of mass momentum. This 
leads to the fact that the FFLO state exhibits a spatially 
varying order parameter. 

In solid state systems, there has been progress toward 
finding experimental evidence of the FFLO state in heavy 
fermion systems [14-.- 17] and also in organic superconduc- 
tors [H, In the context of ultracold gases in one- 
dimensional confinement, there have been experiments 
[8] in qualitative agreement with theoretical studies on 
the FFLO state. However, the question about the exis- 
tence of the FFLO state still remains undecided. This 
subject has received considerable theoretical attention in 
the field of ultracold gases and several experimental pro- 
cedures to probe the FFLO state have been suggested 
to complement the direct imaging of the density profile. 
For instance, the radio frequency (RF) spectroscopy of 
the FFLO state has been a subject of inquiry [20, 21,]. In 
the case of ID systems studies have been made on collec- 
tive mode properties [22], double occupation modulation 
spectroscopy [23] and Josephson junction analogies [24] 
as well as RF specroscopy [25|. Recently, Bragg scat- 
tering and RF spectroscopy were proposed for observ- 
ing the FFLO state in quasi-lD systems [26,]. Moreover, 
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noise correlations have been shown to contain informa- 
tion about the FFLO pairing both in ID and in higher 
dimensions |27l.l28j. 

In this paper we study the density response and collec- 
tive modes of the FFLO state. While several collective 
mode studies exist on imbalanced Fermi gases, only few of 
them consider explicitly the FFLO state [22*, "29", ^SO^ . We 
concentrate on two-compononent spin-imbalanced Fermi 
gases at finite tempereture, in the lowest band of a 2D 
or a quasi- ID optical lattice and with an on-site inter- 
action. The optical lattice is motivated by theoretical 
studies indicating that the lattice aids the formation of 
the FFLO state as the lattice dispersion improves the 
overlap between the Fermi surfaces of the majority and 
minority components [3l|, [32| . Our method for calculat- 
ing the collective mode spectrum is based on the gen- 
eralised random phase approximation (GRPA) for the 
linear response function of the system. The RPA [33| 
is a standard tool for describing collective modes of in- 
teracting fermion systems and it was first applied to the 
BCS context by Anderson [34]. The method has been 
used to describe also e.g layered superconductors [35] 
and the BEC-BCS crossover [3j|, l37|. In addition to 
analysing the collective mode dispersion and the speed 
of sound, we also study the damping properties of the 
collective modes. We find the interesting result that the 
the anisotropic pairing of the FFLO state leads to an 
anisotropy in the speed of sound in the system. Further- 
more, we study a quasi- ID optical lattice in which the 
tunneling in two directions of a 3D lattice is restricted, 
as it has been recently suggested that a quasi- ID geom- 
etry would provide optimal conditions for the formation 
of the FFLO state [38-42]. 

This paper continues in the next section with an in- 
troduction of the Hubbard model and Green's function 
formalism as well as a rederivation of the FFLO Green's 
function. Section ITlI Al outlines the linear response prob- 
lem and the Kadanoff-Baym method [43|, I3 con- 
structing self-consistent linear response approximations. 
We derive the linear response function for the FFLO 
state in section UlIBI After this we present our main re- 
sults in section IIV Al in which we consider the collective 
mode spectrum and the speed of sound in two dimen- 
sional square optical lattices. In section ITVB I we discuss 
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a quasi- ID geometry. Finally, we conclude our work in 
section |Vl 



II. THE FFLO STATE IN A LATTICE 



complex-time /S-matrix is defined as 



S{T,T')=Trexp - / dT"H^{T") 



(4) 



A. The theoretical framework 

We consider a two-component Fermi gas confined to 
the lowest energy band of a square (2D) or cubic (3D) 
lattice with Nl sites, and describe the system with the 
Hubbard model with an on-site interaction. The Hamil- 
tonian Hq of this system is 

ffo = - J.(^.(ri)^i(r2)+Va(r2)^i(ri)) 

(ri,r2),cr 

- ^/i^V^i(r)V^^(r) 

r,cr 

+ Ef/i2V'I(r)V'l(r)^/'2(r)^/'i(r). (1) 

r 

Here a G {1,2} labels the two atomic species, e.q. two 
hyperfine states of a fermionic atom, and are the 
fermionic annihilation and creation operators. (The no- 
tation '0j-(r) is slightly more convenient for the Green's 
function formalism as opposed to the more conventional 
notation Q,cr-) The position vector of the lattice sites is 
denoted by r and the summations run over the set of all 
lattice sites with (ri,r2) meaning summation over near- 
est neighboring sites. Moreover, J is the nearest neigh- 
bour hopping energy, /i^r is the chemical potential and 
Ui2 is the interaction strength between the two species. 
A detailed exposition of the connection of the Hubbard 
model parameters with experimtal parameters for ultra- 
cold gases can be found e.g. in [45]. We employ a periodic 
boundary condition and take the convention h= 1. 

With the assumption that the system is excited by an 
external perturbation i^^ of the form 



The single particle Green's function is then defined as 



E 



,{Ti,T2,t)lpl{Yi)lp^,{Y2) 



(2) 



the total Hamiltonian is H = i^o + i^</). The unperturbed 
Hamiltonian Hq is assumed time-independent, but the 
perturbation i^^ may have an explicit time dependence. 

In the following theoretical treatment we rely on 
Green's function techniques in the Matsubara formalism 
[46, 47] I.e. taking time as a complex parameter, which 
allows us to deal with the finite temperature more effi- 
ciently. The thermodynamic average of the operator O in 
interaction picture in the Matsubara formalism is defined 
as 

0) = — ^ , (3) 



Tr (e-/^^o5(o,/3) 



Here (3 = with the Boltzmann constant and T 
the temperature. is the time ordering operator. The 



G(l,l') = -(T(^(l)V't(l'))). 



(5) 



The shorthand notation 1 is used for the variables ririai. 

It is convenient to extend the range of the spin in- 
dex (J G {1,2} by defining that for a G {3,4} one takes 
ipa- = ~ With this extcusion the 

definition for the Green's function above covers also the 
so called anomalous correlators in which two creation 
or two annihilation operators appear. These functions 
are essential in describing pairing correlations in the sys- 
tem on the mean field level. A similar extension is use- 
ful for J, /i, (j) and U; for ai^a2 G {3,4} one defines 

Jai = ~Jai-2^ Mcri = ~Mcri-2, ^cri,cr2 — ~^cr2 -2,cri -2 ^ud 

Ucri,a2 = U(ji-2,a2-2- thcsc definitions the choice of 
sign allows to write the equations of motion in the most 
fluent form. 

In certain expressions involving two or more field oper- 
ators evaluated at the same time r, the notations r+ and 
r~ specify the time ordering. These notations imply tak- 
ing the limits where r+ ^ r from the positive imaginary 
axis and r~ ^ r from the negative imaginary axis. 

The single particle Green's function follows the equa- 
tion of motion 

5{1, 1') + f (/.(1, 1)G(1, 1') + / E(l, 1)G(1, 1'). (6) 



Here the integral sign is a shorthand notation for sum- 
mation over position and spin in addition to integration 
over time. The overbar indicates a variable of summation 
and integration. Here, the potential (j) appears formally 
as non-local in time but only local potentials are required 
in the work at hand. The inverse non-interacting single 
particle Green's function is 



Go '(1,1') 



+ (7) 



where K^j is the kinetic energy operator defined by 
Kafir) = Ja E(r,r')(/(r') " fi^))- Furthermore, 6(1, 1') 
stands for the Dirac and Kronecker delta for continuous 
and discrete variables, respectively. For a gereneral two- 
body interaction the self-energy S is defined as 

S(l,l') = ni,l)G2(l,i-,2,l+)G-i(2,l'), (8) 

in which the two particle Green's function G2 is given by 
G2(l,2,3,4) = {Tr (^/;(1)V(2)V^(4)V^(3))). (9) 
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For the on-site interaction V{1^ 2) = ^/cricr2^('™i: '™2)^('?"i — 
the self-energy simphfies somewhat due to two trivial 
summations. The equation of motion can also be inverted 
for as 

1') = Go\l, 1') - HI, 1') - !')• (10) 



B. Imbalanced superfluid in a mean field model 

In this section we consider the FFLO mean field de- 
scription for a spin- imbalanced Fermionic superfiuid first 
studied by [ll|, |l2j . Notice that in this section the focus 
is on ground state properties and the external field i^^ is 
not needed. 

The equation of motion for the single particle Green's 
function involves the two particle Green's function (|9]). 
This object follows again its own equation of motion 
which involves the three particle Green's function, and 
continuing this way, one derives an infinite set of equa- 
tions known as the Martin- Schwinger hierarchy. In prac- 
tise, one needs to decouple this hiearchy on some level in 
order to find the single particle Green's function. We re- 
sort to a mean field approximation in equation (j6j). This 
is often referred to as the Hartree-Fock-Gor'kov approx- 
imation. In this approximation and with the notation 
U = Ui2 the self-energy is 



"— G22 


G12 





— G14" 


G21 


-Gil 


—G23 








— Gs2 


— G44 


G34 


_— G41 





G43 


— G33_ 



• (11) 



All of the Green's functions appearing in S have the vari- 
ables (riTi, riT;^). 

In a system with uniform density, the Hartree terms 
on the diagonal can be absorbed into the chemical po- 
tentials. The Fock-exchange terms with G12, G21, G34 
and G43 are negligible in solving for the ground state 
for ultracold Fermi gases, as spin-fiips are energetically 
highly unfavourable in experimentally relevant magnetic 
fields. Finally we introduce the key element of the mean 
field FFLO theory. We assume that the pairing correla- 
tions have an oscillating structure so that the self-energy 
is 



^=S{riri,r[r[)A 







-2iqTi 






-2iq-ri 







o2iq-ri 






o2iq-ri 







(12) 



Here q is the FFLO pairing vector. In the special case of 
q = and Ni = N2 the system is reduced to the standard 
BCS description. The case with q = and A^i 7^ N2 is 
commonly known as the breached pair (BP) state. 

The quantity Ae"^*^ '^^ is the order parameter of the 
FFLO state. In general, A is related to the energy gap of 
pair breaking excitations. We point out that this choice 



of the order parameter is not the only possibility, and it 
has been shown theoretically [1^ [l3| that for instance 
an order parameter of a cosine form would be energeti- 
cally more favourable. However, the current choice allows 
for developing the theory analytically much further, thus 
making the physics more transparent. 

To quarantee the consistency of the mean field solu- 
tion, we must have 



f?32(l,l+) = ^Ae-2*i-\ 



(13) 



This condition is the FFLO gap equation. In the FFLO 
self-energy all the nonzero elements are connected by 
complex conjugation or anticommutation relations of the 
field operators. Therefore, one indeed has just one inde- 
pendent gap equation. 

One can also fix the expected particle numbers for each 
atom species, A/'^r, with the number equations 



iV. = ^G,.(fr,fr+). 



(14) 



However, for a uniform density distribution one may 
write the number equation directly in terms of the den- 
sity (or more precisely the filling fraction) Ua- = N^j/Nl. 
The number equation is 



V = Gaa(rr,rr+). 



(15) 



In order to find out the values of A, jii and 112 for any 
given FFLO state with pairing vector q and particle num- 
bers Ni and N2 we need to find a solution for the gap 
and number equations for the state. 

In the following we derive a closed algebraic form for 
the FFLO Green's function in momentum and frequency 
space and rewrite the gap and number equations accord- 
ingly. 

It is possible to solve the Green's function in the 
present approximation analytically in the Fourier space. 
Here, in order to find an algebraically closed set of 
Fourier components we have to pay particular atten- 
tion to the broken translation invariance of the FFLO 
order parameter. However, our system is still trans- 
lation invariant with respect to time. Thus, we have 
G(ri, r2) = G(ri —T2) and we may take the Fourier trans- 
formation in time directly with respect to ri — r2 . 

We define the Fourier transformation of the Green's 
function as 



G(pi,P2,^) = ^ j d{Ti-T2)i 



,iuj{Tl-T2) 



T{pi • ri)G(ri,r2,ri - r2)T\p2 • 1*2), 

(16) 



where the Fourier transfrom matrix is given by 



^(Pi -ri) 



g-*Piri 

e-^Pi'^i 

e^Pi'^i 

e^Pi-^i 



. (17) 
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Here pi and p2 are momenta and is a frequency (or en- 
ergy, as we have chosen the convention h = 1). The sign 
convention of has been chosen so that it agrees with 
the Fourier transformation of the field operators. Due 
to the periodic boundary condition in complex time, also 
the frequency spectrum is discrete covering the fermionic 



Matsubara frequencies uj 



(2n+l)7r 



where n is an integer. 



The inverse Fourier transformation is 
1 



G(ri,r2,ri - r2) 



E 



e-^"^"^-"^^^^(pi • ri)G(pi,p2,c^)^(p2 • r2). (18) 

Here the momentum summations run over the discrete 
momentum spectrum and the frequency summation over 
the Matsubara frequencies. 

We now Fourier transform the equation of motion ([6]) 
for G in the FFLO state. For brevity, we deal first with 
the (Ji, (72 G {1, 4} block, which has the Fourier transform 



iuj - ^i{pi) 

za; + ^2(pi) 



^Pl,P2^ 



1 

1 



G(pi,p2,a;) 
G(2q-pi,p2,c^). (19) 



Here the non-interacting particle energy ^cr(p) is 

^^(p) = e(p)-/i,, (20) 

in which e(p) is the lattice dispersion. 

Equation ([19]) is mixing different Fourier components 
of the Green's function due to the broken translation in- 
var iance of the order parameter. To be more explicit, for 
Gil we have 



(icj -^i(pi))Gii(pi,p2,cj) 
=^Pi,P2 + AG4i(2q - pi, p2,a;). 



(21) 



However, the equation is still closed. Relabelling pi with 
2q — Pi the equation for G41 (2q — pi , p2 , Co') is 

{iu + ^2(2q - pi))G4i(2q - pi, p2, cj) 
=AGii(pi,p2,a;). (22) 

From this pair of equations it is straightforward to solve 
for Gil and G41. In a similar manner we find the Green's 
functions G14 and G44. The solution is 



G'ii(pi,P2,cj) Gi4(pi,2q- P2,cj) 

G4i(2q- pi,p2,a;) G44(2q - pi, 2q - p2, a;)_ 

^Pl,P2 

" (iw-a(pi))(*w + e2(2q-pi)) - A2 

iw + 6(2q-pi) A 

A iw-^i(pi) 



(23) 



o'i,cr2 G {2, 3} is similar and can be written as 



G22(2q - Pi, 2q - p2,cj) G23(2q - pi, P2, cj) 

G32(pi,2q- P2,CJ) G33(P1,P2,CJ) 

^Pl,P2 

" {ioj + a(Pi))(iw - e2(2q - pi)) - A2 

iw + Ci(pi) -A 

-A iw-^2(2q-pi) 



(24) 



Notice that all the other Green's functions Gau are triv- 
ially zero in the adopted approximation. 

The solution can be written in a form which is easier 
to analyse and is similar to the conventional form for 
the BCS Green's functions. This final step makes the 
application of finite temperature Matsubara summation 
techniques straightforward. One defines the well-known 
quasiparticle energies E± pT] as 



E±{p) = ± 



6(p)-6(2q-p) 



a(p)+6(2q-p) 



A2, 



and the coherence factors u and v as 
u{p) = 

u(p) = 



/ E+(p) + e2(2q-p) 
E+{p) + E_ip) ■' 

/i?_(p)-6(2q-p) 



E+{p)+E_{p) 



(25) 

(26) 
(27) 



With these definitions the FFLO Green's functions can 
be written for the block cri,(j2 G {1,4} as 



G'ii(pi,P2,^) Gi4(pi,2q- P2,CJ) 

G4i(2q- pi,p2,cj) G44(2q- pi,2q- p2,a;) 



^Pl,P2 



iuj - E^{pi) 



^(Pl)^ ^(Pl)^(Pl) 
^(Pl)^(Pl) ^(Pl)^ 

V(pi)2 -^(Pl)^(Pl) 
-^(Pl)^(Pl) ^(Pl)^ 



iuj -I- ^-(pi) 
and for the block ai, (J2 G {2, 3} as 



(28) 



G22(2q- pi,2q- p2,cj) G23(2q - Pi, P2, cj) 
G32 (pi , 2q - p2 , cj) G33 (pi , P2 , cj) 



^Pl,P2 



iuj - £^-(pi) 

^Pl,P2 



u{pif -u{pi)v{pi) 

-^(Pl)^(pl) ^(Pl)^ 

^(Pl)^ ^(Pl)^(pl) 
^(Pl)^(pl) ^(Pl)^ 



(29) 



The solution for the Green's functions in the matrix block 



While the normal Green's functions are diagonal in mo- 
mentum space, the anomalous Green's functions are not, 
reflecting the oscillatory structure of the pairing field. 

Finally, we present the gap equation (p!3|) and the num- 
ber equations (p!^ in Fourier space. The inverse Fourier 
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transform of G32(l, 1^) appearing in the gap equation is 
G32(rr,rr+) = 

^ e-'^^-u{pMp){l - n^E^ip)) - n^(^_(p))), 

(30) 



Nl 



in which the Fermi distribution np has been obtained 
from the Matsubara summation 



1 ^ giu;(r + -r) 



(31) 



Thus, the gap equation is 
A _ 1 



u{p)v{p){l - nF{E+{p)) - (i5_(p))). 

(32) 



Similarly the number equations in terms of the filling 
fraction are 

ni = Gii(rr,rT+) = 
^-Y,HpfnF{E+{p)) + v{p)Hl - np{E.{p))), 



N, 



(33) 



and 



n2 G22(rr,rr+) 
^-^u{pfnF{E_{p)) + ^(p)^(l - n^(^+(p))). 



Nl 



(34) 



III. DENSITY RESPONSE IN THE FFLO 
STATE 

A. Linear response theory 

Having estabhshed a ground state description for the 
imbalanced superfluid, let us turn to study density fluc- 
tuations caused by potentials that couple to the particle 
density. The Hamiltonian i^^ for such external potentials 
is of the form 

+ 1 Mrt)4{^)Mr)- (35) 

Now, if and when solving the system with H = Hq + H^p 
directly is not feasible, progress can be made by assuming 
that Hff) is a small perturbation. To find out the effect of 
H(h for instance on the density of atom species a = 1 one 



can then write ricr-^ (rt) as a variational series with respect 
to (j) and truncate this series to the first order, obtainging 



ni 



(r.) = (n.(r,),.„./^.(.-)(|||;^^^ 



(36) 



Let us continue in the Matsubara formalism. The con- 
cept above can be generalised for any Green's function 
by writing it as a variational series to the first order with 
respect to </)(!, 2) so that 

SG{1,1'] 



G(l,l') = (G(l,l'))^^o+ / 0(2,3) 



<5^(2,3) /^=o 

(37) 



Notice again, that the overbar indicates summa- 
tion/integration over position, time and spin. The varia- 
tional derivative in the equation above defines the linear 
response function 



^(''''')-lM2^2) 



(38) 



= 



For example the density is given by ricr^ (riri) = G(l, 1+) 
so L{12^ 1+2') would give the density response function. 

The linear response function carries information about 
the excited states of the unperturbed system. If one is 
able to solve the linear response function one can then 
extract from it the excitation spectrum of the system. 
For instance the collective density modes appear as sim- 
ple poles of the density response function in frequency 
space. 

Following Kadanoff and Baym [isl, |44] one can derive 
an equation for the linear response function from the 
equation of motion (p!Q|) . To outline the derivation of 
their result briefiy, let us begin by taking the variational 
derivative of the identity J GG~^ = 6^ which yields 



SG{1,1') 
(50(2', 2) 



Now, inserting G ^ from (p!Q|) into the previous equation, 
evaluating the expression at = and identifying L we 

get 

Since the self-energy does not depend explicitly on the 
external perturbation, the chain rule of differentiation 
gives 

<5E(3,4) f ,5S(3,4) 5G(5,6) 



54>{2',2) 



SG{5,6) (50(2', 2) 
(5E(3,4) 



(5G(5,6) 



L(52,62'). 



(41) 
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Therefore, one may write equation (|4Q|) explicitly as an 
integral equation for the linear response function 

L(12,1'2') = 

+ /G(l,3),.oG(4,lO,.o(g||)^_^i(52,62'), 

(42) 

which is the result of [isl, |4^. This equation for the 
response function quarantees a self-consistent theory in 
the sense that the linear response function obeys the same 
conservation laws as does the single particle Green's func- 
tion. From this point on, we leave out the notation ^ = 
as in the following all of the variational derivatives are 
evaluated at (/) = 0. Notice that in the equation above 
there is a trivial variational derivative 



<^0(3,4) 
50(1,2) 



5{Tzn - riTi)(5(r4T4 - r2T2) (5 ((73(74,(71(72), 



(43) 



where we have the notation 



J ((73 (74, (71(72) = 

(5((73,(7i)(^((74,(72) -(^((74,(71 ±2)(^((73,(72 ±2). (44) 

The second term in this definition owes to the fact that 
for the extended index (7i,(72 G {3,4} we have 0^1, cr2 = 

— (/>cr2-2,cri-2. 



B. Derivation of the FFLO density response 
function 



We derive the FFLO density response function in this 
section. The Kadanoff-Baym method applied to the 



Hartree-Fock-Gor'kov approximation (pTj) leads to the 
generalised random phase approximation (GRPA). We 
now extend this to the case of a FFLO state. In the 
density response problem, several simplifications to the 
general equation for the linear response function are ap- 
parent. First of all, it is sufficient to consider a local per- 
turbation of the form (/)(!, 2) = (/>cricr2 ('^i'^i)^('^i'^i ~'^2'?"2). 
Furthermore, the density operator itself as well as all the 
Green's functions in the self-energy are local. Thus, the 
density response function of which we are interested in 
contains only terms with ri = r^, ri = t[ and Y2 = 1*2? 



r2 



Moreover, due to the time translation invar iance 



this response function depends only on the time differ- 
ence Ti — r2 . We then find it suitable for our purposes to 
define the notation 

^aia2a;a^(ri,r2,ri -T2) = 

(5(riri, r;rO(5(r2r2, 4^^)^(12, lY). (45) 

The density response functions |^ and |^ in equation 
(|36l) correspond to Lim and L1212, respectively. 
Equation (|4Q]) now reads 

^aia2a;a^(ri,r2,ri - r2) = 

cr3iri,r2,Tl T2) 

X G'^4a;(r2,ri,r2 - n)^ (0-3 0-4, 0-3(72) 
Ga,^3(ri,f3,ri - fs) 

X G'^4a;(r3,ri,f3 - ri)^^^^(r3, r2, r3 -r2). (46) 



Let us now write down explicitly the equation for Lim 
to show how to proceed with the solution. Since many 
of the FFLO Green's functions are identically zero, just 
one term remains from the direct coupling to the exter- 
nal perturbation which is the first term in equation (|^6]) . 
Similarly, only four non-zero terms arise from the second 
term in equation (|^ . 



Liiii(ri,r2,ri - r2) = Gii(ri,r2,ri - r2)Gii (r2, ri, T2 - n) 





1 Gi4(ri, 


r3, 


n 


- ^3)^41(1^3, 




rs 


- ri)Liiii(f3, 


1*2, 


T3 


-T2) 


"J 


f Gii(ri, 


r3, 




- r3)G4i(r3, 


ri, 




- ri)Lii4i(f3, 


1*2, 


T3 


-T2) 


"J 


f Gi4(ri, 


^3, 


n 


- r3)Gii(r3, 






- ri)^4iii(f3, 


1*2, 




-r2) 


u i 


f Gii(ri, 


r3, 




- T3)Gii(f3, 






- Ti)^414l(r3, 


1*2, 




-r2) 



I 

In forming this equation we have used the identities ^^3131 = —^1111 and L2121 = — ^4i4i- This substitution 



7 



owes to a more general identity: Noting the definitions 
(|38|) of the response function and (|5]) of the Green's func- 
tion one concludes based on equal time anticommutation 
relations that 



^cr,a,zy,/3(ri,r2,r) = — I/^+2,a,zy+2,/3(ri, r2, r), 

^a+2,a,z.,/3(ri,r2,r) = -i^a,a,i.+2,/3(ri,r2,r). (48) 



As with the Green's function, the integral equation 
for Lull (|47|) can be cast into an algebraic equation 
in Fourier space. Notice that for instance in spherically 
symmetric harmonic trapping geometries one can make 
similar progress with the choice of the harmonic oscillator 
states as basis functions [48]. We use the same Fourier 
transformation convention for Lcr,i,z/,i = we de- 

fined for Gfju in ([11]) ^-e. the sign convention is given by 
a and u. The Fourier transformation yields for Lim, i.e. 
the equation 



^iiii(pi,P2,^) = (^pi,p2niiii(pi,a;) 
+/7ni44i(pi,cj)Liiii(pi,p2,cj) 
-UUii4i (pi , a;)Lii4i (2q + pi , -p2, cj) 
-[/ni4ii (pi , a;)L4iii (2q - pi , p2 , cj) 
-hC/niiii(pi,a;)L4i4i(-pi, -P2,cj). 



(49) 



Here we have the notation 

n 



0-10-20-30-4 (P^ ^) 

^^1^2 (V (P + s), (P + s), X + cj) 



xGa3a4(Aa3(s), (s),x), 

where Acr(p) is defined so that 



A.(p) 
A.(p) 



P, 
2q 



P, 



'^£{1,2}, 
<JG{3,4}. 



(50) 



(51) 



Equation (|47|) is rather analogous to the previously solved 
equation for the FFLO Green's function. We see that 
we need to construct equations also for L1141, L4111 and 
L4141 . Using again the symmetry property (j48j) for these 
equations we find out that no other linear response func- 
tions enter the equation. The final task is to identify 
those Fourier components which form a closed equation. 
With this rationale one arrives at the following matrix 
equation 

M(i)(pi)L(i)(pi,p2) = ,5p„p,n(i)(pi). (52) 

Here the vector of linear response functions L^^^ is defined 
as 



(53) 



^llll(Pl,P2) 

^Ii4i(2q + Pi,-P2) 
L4iii(2q- pi,p2) 

L ^414l(-pi, -P2) J 

On the right hand side 11^^^ contains the terms from the 
direct coupling to the external perturbation 



n(i)(pi) 

The coefficient matrix M^^^ is given by 



■niiii(pi)" 
niii4(pi) 
n4iii(pi) 
n4ii4(pi). 



(54) 



J 



M^i)( 



PiJ 



-ni44i(pi) nii4i(pi) ni4ii(pi) -niiii(pi)" 

-ni444(pi) nii44(pi) ni4i4(pi) -niii4(pi) 

-n444i(pi) n4i4i(pi) n44ii(pi) -n4iii(pi) 

-n4444(pi) n4i44(pi) n44i4(pi) -n4ii4(pi). 



(55) 



One derives similarly for L1212, i-e. 4^, the equation and n*^^^ is 



M(2)(pi)i:(2)(p,,p2) 

in which L^^^ stands for 



^Pl,P2 



n^'Hpi), (56) 



i^^'^(pi,P2) 



^1212(Pl,P2) 

^I242(2q + Pi,-P2) 

^4212(2q - Pl,P2) 
. ^4242( — Pi, — P2) . 



(57) 



n(2)(pi) = 



■ni44i(pi)' 

ni444(Pl) 

114441 (Pl) 

n4444(pi). 



(58) 



The coefficient matrix M^^^ turns out to be the same as 
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In our theoretical treatment we have an obvious sym- 
metry with respect to interchanging indices 1 and 2 and 
indices 3 and 4. Therefore, we obtain the equation for 
the density response functions L2222 and L2121 directly 
from the results above. 

The frequency summation in the definition of 11, equa- 
tion (|5Q|) , can be handled analytically. One applies the 
identity 



-y 



1 



1 



El 



(59) 



to the four cross terms that arise when one inserts the 
Fourier transformed Green's functions ([28|) and ([29|) to 
equation (|5Q]) . The momentum summation in (|5Q|) needs 
to be carried out numerically after which one simply in- 
verts the matrix equations presented above to obtain the 
Matsubara, i.e. imaginary frequency, response function. 
From this, the real frequency retarded linear response 
function is obtained by means of analytical continuation. 



IV. RESULTS 
A. 2D square lattice 

In the following results we consider a system in a two 
dimensional square optical lattice with lattice constant d 
and Nl N^Ny 40000 lattice sites with N^j, = Ny = 
200 lattice sites in each direction. We assume that the 
perturbation potential is the same for both atom species, 
i.e. (pi = (j)2 in which case it is most natural to study 
the density response function xi(k, co') = Liiii(k, cj) + 
^1212 (k,a;) where Lmi and L1212 are the responses of 
density ni to potentials (j)i and 02, respectively. In the 
following, the parameters are chosen so that ni is the 
density of the majority component. Similar conclusions 
hold also for the response of the minority component. 
The assumption of (j)i = (j)2 is not a crucial one as the 
collective mode dispersion is the same for all choices of 
these potentials. Moreover, we assume that the FFLO 
vector q is directed along the x-axis. 

We solve equations ([52]) and (|56l) numerically for imag- 
inary frequencies. We then use a Fade approximant [49] 
to carry out the analytical continuation and obtain the 
response function for complex frequencies. 

In Figure [T] we plot the real part of the density re- 
sponse function for typical parameters as a function of 
the wave vector k and the real frequency uo. In this fig- 
ure we have chosen k parallel to q. The collective density 
mode appears as nearly diverging feature of the density 
response at particular a wave vector k and frequency uj 
Towards higher wave vectors we see a typical broadening 
of this feature owing to the increase of damping. The 
actual eigenfrequency of the mode is complex with the 
frequency u — ij where 7 is the damping rate. For the 
small damping rates the response has a very sharp jump 




k (units of 1/d) 



FIG. 1. (Color online) The real part of the density response 
Xi(k,c<;) for typical parameters /ii = 3. 5 J, /i2 = 2. 5 J, T = 
0.07J, U = -3.0J, A = 0.27J and q = 127r/((iiVa.) with k 
parallel to q. The collective mode is seen as a clear divergent 
behaviour of the density response. The lone peak on the kx 
axis is a numerical instability. 



also as a function of real frequencies. The collective mode 
is gapless with a linear dispersion for small k. It corre- 
sponds to the Anderson-BogoHubov phonon of the BCS 
state of a neutral Fermi gas [34]. 

We solve the dispersion relation co'(k) and the damping 
rate 7(k) by solving the poles of the response function xi- 
We plot cj(k) and 7(k) in Figure [2] for wave vectors along 
the X- and y-axes i.e. parallel and perpendicular to the 
FFLO vector q. Notice that the numerical method pro- 
duces several instabilities in the damping rate at higher 
wave vectors while the real frequency uj is far more sta- 
ble. The speed of sound to the direction of the x-axis, 
can be obtained from the dispersion by the definition 



duoikex) 



dk 



(60) 



where e^; is the unit vector in x-direction. One defines 
Cy similarly. For the dispersions presented in Figure [2] 
we obtain Cx = l.39Jd and Cy = 1.2AJd. We observe 
that the finite FFLO vector causes a clear difference be- 
tween the parallel and perpendicular (w.r.t. q) speeds of 
sound. In this case the relative difference in the sound 
propagation is Cx/cy — 1 = 12%. This observation sug- 
gests immediately an experiment in which one creates a 
local density perturbation in the system and then mon- 
itors the propagation of this perturbation to collect in- 
formation about the FFLO state. Such an experiment 
would create a rather remarkable contrast with e.g. a 
breached pairing state, which is isotropic with Cx = Cy. 

Turning back to analyse the damping rates presented 
in Figure [21 we find that the damping rate is only a 



9 



1 

0.9 

0.8 

0.7 

2 0-6 
o 

0.5 
t 0.4 
0.3 
0.2 
0.1 



ooo 



oo 



oooOq" Oo 



❖ y 



0.2 0.4 0.6 0.8 1 

k (units of 1/d) 



1.2 





1 




0.9 




0.8 




0.7 




0.6 


"o 




w 


0.5 


"c 








3 


0.4 




0.3 




0.2 




0.1 








, i 11 !! I 



0.2 



0.4 0.6 0.8 

k (units of 1/d) 



1.2 



0.16 
0.14 
0.12 

o 

^ 0.08 
c 

^ 0.06 
0.04 
0.02 



yy 



^ 0^. ^ ^ A ^ ^ <^ 
^A <^ A ^ 

^ O 

g doooo 
o 



oO 



0.2 0.4 0.6 0.8 

k (units of 1/d) 



1.2 



FIG. 2. (Color online) The dispersion relation C(;(k) and the 
damping rate 7(k) calculated for /ii = 3.5 J, /J2 = 2.5 J, T = 
0.07J, U = 3.0J, A = 0.27J and q = 127r/((iAra:) along the 
X-axis (circles) and y-d^xis (diamonds). The FFLO vector q is 
directed along the x-axis and this anisotropy creates a clear 
difference between the two dispersions and damping rates. 



fraction of the mode frequency, with jx/^x = 5% and 
Jy/ojy = 7% for the lowest frequencies. With this ob- 
servation we conclude that the collective modes can be 
considered as well-defined elementary excitations of the 
FFLO state. An important aspect here is that the damp- 
ing is fundamentally not an effect caused by the finite 
temperature. In contrast to the BCS state, the FFLO 
state contains also at zero temperature unpaired quasi- 
particles which cause a finite lifitime for the collective 
modes. In Figure [3] we illustrate the locations of the 
quasiparticle channels for the system of Figures [1] and [2j 
Indeed, the collective mode dispersion lies deep within 
the region containing quasiparticle transitions. To com- 
pare, the pair breaking excitations are present also in a 
BCS state of a neutral Fermi gas, but at zero tempera- 
ture the quasiparticle transitions are completely absent 



FIG. 3. (Color online) An illustration of the quasiparti- 
cle transitions for a system with /ii = 3. 5 J, /i2 = 2. 5 J, 
A = 0.27J and q — 12Ti/{dNx) along the x-axis, correspond- 
ing to a zero temperature Fermi distribution. A cross denotes 
a transition of a single quasiparticle and a plus a pair break- 
ing which involves the creation of two quasiparticles. Notice 
that the figure shows only the positions of the quasiparticle 
excitations, not their relative weights. The small void regions 
along the kx axis owe to the finite system size. 



and the Anderson-Bogohubov phonon is thus undamped. 

In Figure [4] we study the speed of sound as a function 
of the length of the FFLO vector q when the other sys- 
tem parameters are held constant and the gap is solved 
from the gap equation (p!3|) for each q. In the range 
qdNx/ (27r) = 4 ... 8 we find an anisotropy of 5 to 7 per- 
cent. The figure also indicates that the speed of sound 
tends to increase with q. This is caused by the fact that 
the gap A decreases when q increases, which in turn 
causes the increase in the speed of sound. The two small- 
est possible FFLO vectors have been left out, since they 
lead to a pair breaking gap which is larger than the chem- 
ical potential difference 2 A > \fii — fi2\ and therefore the 
state is not physically relevant. 

We then turn to examine the effect of the polarisation 
P — {N1—N2) I {Ni-\-N2) on the speed of sound in Figure 
[5l Technically, we vary the chemical potential difference 
while holding the average chemical potential (/ii + /i2)/2 
and other system parameters constant. We find that the 
anisotropy in the sound propagation increases from 2% to 
10% when the polarisation increases to the value P = 0.2. 
Also the speed of sound itself increases with the polarisa- 
tion, and the explanation is analogous to the discussion of 
the FFLO vector above: while holding the other system 
parameters fixed, the gap A decreases with the chemi- 
cal potential difference, or the polarisation. In the range 
P < 0.11 the pair breaking gap is larger than the chem- 
ical potential difference, and the FFLO ansatz does not 
produce a physically relevant state. 

We examine the tempereture dependence of the speed 
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FIG. 4. (Color online) The speed of sound parallel (cx, circles) 
and perpendicular (cy, diamonds) to the FFLO vector q = 
QBx as a function of length of the FFLO vector q for a system 
with fii = 3.4J, fi2 = 2.6J, T = O.IJ, U = -2.8J. The gap 
A is solved from the gap equation (p^ for each q. 




FIG. 5. (Color online) The speed of sound parallel (cx, circles) 
and perpendicular (cy, diamonds) to the FFLO vector q = 
qGx cLS a function of the polarisation P {Ni-N2)/{Ni-\-N2). 
Here we vary the chemical potential difference S/j^ — fii — /J2 
to vary the polarisation while holding (/xi + /i2)/2 constant 
at 3.0 J. The other parameters for the system are T = O.IJ, 
U = -2.8J, q = 127r/{dNx). The gap A is solved from the 
gap equation ()13p for each set of parameters. 



of sound in Figure [6l Firstly, we stress that while the 
temperature does affect the dispersion relation, it always 
remains true that in the small wavelength limit the col- 
lective mode is massless and the dispersion goes linearly 
to zero. We see that the speed of sound and in particular 
the anisotropy is fairly robust against changes in tem- 
perature deep in the superfluid phase. When the tem- 
perature approaches the critical temperature (in Figure 
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FIG. 6. (Color online) The speed of sound parallel (cx, circles) 
and perpendicular (cy, diamonds) to the FFLO vector q = 
qex as a function of the temperature T. Here the system 
parameters are in the figure above /xi = 3.4J, /J2 = 2.6J, 
U = — 2.8J and q — 127v/{dNx) and for the figure below 
fii = 3.5J, fi2 = 2.5J, U = -3.0J and q = 127i/{dNx). Note 
that the critical temperature when holding these parameters 
constant is Tc = 0.15 J (above) and Tc = 0.12 J (below). 



[6]Tc = 0.15J (above) and = 0.12 J (below)) the speed 
of sound increases. This effect arises from the tempera- 
ture dependence of the gap A which falls to zero when 

In Figure [71 we plot the damping rates of the k = 
27T/{dNx) collective mode. This is the lowest non-zero 
wave vector. As mentioned previously, in the FFLO state 
finite damping is present already at T = in contrast to 
a BCS state, since the FFLO state has quasiparticle exci- 
tations even at T = 0. Increasing the temperature creates 
more quasiparticle excitations and the damping becomes 
stronger. The result implies that in a broad range of 
temperatures a significant portion of the damping can 
be attributed to the zero temperature quasiparticles as 
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FIG. 7. (Color online) The damping rate 7 of the lowest k- 
mode, k = 27r /(dNx) in directions parallel (7^, circles) and 
perpendicular (7^, diamonds) to the FFLO vector q = qex 
as a function of temperature. The parameters for the plots 
above and below are the same as in figure (6] 



opposed to thermal excitations. The damping is notably 
enhanced close to Tc as the excitation gap A decreases 
rapidly close to Tc. The largest damping rates in the 
Figure [3 are 13% (above) and 15% (below) of the cor- 
responding collective mode frequency, which means that 
the low frequency collective modes are well-defined in the 
entire temperature range shown. 
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k (units of 1/d) 



FIG. 8. (Color online) The real part of the density response 
function Liiii(k, uj) for an FFLO gas in a quasi-lD optical 
lattice along the x-axis. The hopping in the y and z directions 
is Jy = Jz = 0.02Jx and moreover /ii = 2.6Jx, /i2 = lAJx, 
T = 0.05 Jo:, U = -2.7 Jx, A = 0.3 Jx and q = Stv/ (dNx) 
with q parallel to the x-axis. There are two strongly peaked 
responses. The lower branch corresponds to the phonon mode 
of the previous section while the upper branch is caused by 
quasiparticle excitations centered into a narrow stripe due to 
the quasi- ID nature of the system. 
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B. Quasi-ID lattice 

We now turn to quasi- ID optical lattices motivated by 
recent theoretical works which predict that the FFLO 
state would be more stable in such geometries [38'-'42j. 
We study a three dimensional cubic lattice with lattice 
constant d and a weak hopping in two directions i. e. Jy = 
Jz ^ Jx- The FFLO q- vector is directed along the x- 



FIG. 9. (Color online) The zero temperature quasiparticle 
channels (for quasiparticles with zero transverse momentum) 
calculated for the system of Figure [51 A cross denotes a transi- 
tion of a single quasiparticle and a plus a pair breaking which 
involves the creation of two quasiparticles. At low wave vec- 
tors there is only a narrow stripe of quasiparticle channels 
which becomes clearly visible also in the GRPA response. 
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axis. 

In the quasi- ID setup the collective mode spectrum 
along the different axes scales with the appropriate hop- 
ping strength. Therefore, one observes that Cyjc^ ~ 
JylJx < 1. Now, it is stih true that the FFLO 
state creates an additional difference perpendicular to 
the q-vector but the relative variation in Cy and is 
in the same order of magnitude as in the case of an 
isotropic lattice and therefore negligible in comparison to 
the anisotropy caused directly by the hopping anistropy. 
Therefore, the analysis of the sound propagation in differ- 
ent directions does not suggest direct means for observing 
the FFLO state in a single experiment in systems that 
are strongly anisotropic already by geometry. 

However, there is another interesting feature we find 
at the quasi- ID limit. In the following we have an op- 
tical lattice with Nj^ = 120^ lattice sites i.e. Nx = 
Ny = Nz = 120 lattice sites in each direction, and 
Jy = Jz = 0.02Jx. In Figure [8] we plot the density 
response function Liiii(k, cj) for an FFLO state with 
fii = 2.6J^, fi2 = lAJx, T = 0.05J^, U = -2.7J^, 
A = 0.3Jx and q = 87r/{dNx). We now find two strongly 
peaked responses. Both of these branches are linear 
for small wave vectors and also gapless. Moreover, the 
higher branch vanishes for kd ~ 0.5. The lower branch 
corresponds to the collective density mode of a multidi- 
mensional lattice. The higher branch in turn is created 
by the FFLO quasiparticle transitions which are illus- 
trated in Figure [9l This higher branch is present also in 
L1212 but with the opposite sign. Therfore it tends to 
cancel, though not exactly, the contribution of Lmi in 
Xi = ^1111 + ^1212- To illustrate the effect clearly we 
thus plot in this section Lmi. 

In the quasi- ID limit the quasiparticle energies are 
dominated by the kx wave vector and the transverse mo- 
mentum only creates a relative variation on the order of 
Jy/Jx to these energies. This ID-like dispersion of the 
quasiparticles places a notable restriction to the possi- 
ble energy and momentum trasfers for small energies and 
momenta. On the other hand, each possible quasiparticle 
transition gains a significant weight for the very same rea- 
son by the following argument. Let SEo{kx) be the energy 
of transition between an empty and a filled quasiparti- 
cle state with no transverse momentum. This means that 
there are filled and empty quasiparticle states at energies 
E±{kx ^Px.Py = 0,pz = 0) and E±{px,Py = 0,Pz = 0) 
for which dEo{kx) = E±{kx^Pa:, 0,0)- E±{px, 0,0). De- 
noting then SE{kx) = E±{kx^Px,Py,Pz)-E±{px,Py,Pz) 
we find after some straightforward algebra that 5E{kx) — 
dEo{kx) ~ Jykxd, for small kx- In other words all the 
quasiparticles with different momenta Py and Pz but the 
same momentum Px contribute at a narrow energy inter- 
val on the order of Jykxd and therefore we find a strongly 
peaked response in Figure [8] following the quasiparticle 
energies of Figure O The strong quasiparticle response 
vanishes towards higher wave vectors as the quasiparticle 
transitions are spread to a wide energy range. Moreover, 
for low wave vectors the dispersion of the lower collective 



mode does not overlap with the quasiparticle excitations 
and therefore the mode is undamped, unlike in the higher 
dimensional case. 

To compare to the results of [22] for a strict ID sys- 
tem and an LO ansatz (cosine form order parameter), it 
is interesting to notice that such a system has a similar 
two mode structure at low wave lengths as the FF ansatz 
(plane wave order parameter) in a quasi- ID setup. How- 
ever, at higher wave lengths the mode associated with 
the excess quasiparticles does not exhibit an additional 
Brillouin zone structure in the FF case. The reason is 
that this Brillouin zone structure is caused in the LO 
case by the oscillatory structure of the quasiparticle den- 
sity profile [22,] , while the FF quasiparticle density profile 
is uniform. 



V. CONCLUSIONS 

We have studied the density response of a spin- 
imbalanced ultracold Fermi gas in an optical lattice in 
the FFLO state. Using the Kadanoff-Baym formalism we 
derived the linear response function for this system in the 
generalised random phase approximation. We then calcu- 
lated the collective mode spectrum in a 2D square optical 
lattice and showed that the speed of sound is anisotropic 
due to the anisotropy of the FFLO pairing. This suggests 
an experiment in which one monitors the propagation of 
a local density perturbation in order to find evidence of 
the anisotropic pairing mechanism of the FFLO state. 

Moreover, we studied the damping of the collective 
modes and showed that despite the presence of quasi- 
particles in the FFLO ground state the collective modes 
have a relatively weak damping rate and are thus well- 
defined and physically meaningful elementary excitations 
of the system. 

We also studied a quasi- ID system. In this case the 
anisotropy of the sound propagation is predominantly 
caused by the anisotropy of the lattice itself in contrast to 
a possible exotic pairing mechanism. However, the quasi- 
ID system is qualitatively different from higher dimen- 
sional systems as it contains an additional collective-type 
response of quasiparticles. 

To draw future quidelines, a clear way to improve on 
the results presented in this paper would be the inclusion 
of more elaborate order parameter structure, in partic- 
ular the LO ansatz with a cosine type order parameter. 
In this case, one cannot simplify the problem in momen- 
tum space to the same extent as in section UlI Bl of this 
paper, and a heavier numerical method in position space 
is called upon. One would still anticipate an anisotropic 
speed of sound for the LO ansatz as well, based on pre- 
senting the state as a superposition of two FF states, as 
well as by the arguments of [29]. 
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